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ABSTRACT 

We investigate the nature and spatial variations of turbulence in the Small Magellanic Cloud (SMC) 
by applying several statistical methods on the neutral hydrogen (HI) column density image of the SMC 
and a database of isothermal numerical simulations. By using the 3rd and 4th statistical moments we 
derive the spatial distribution of the sonic Mach number (A4 S ) across the SMC. We find that about 
90% of the HI in the SMC is subsonic or transonic. However, edges of the SMC 'bar' have A4 S ~ 4 
and may be tracing shearing or turbulent flows. Using numerical simulations we also investigate how 
the slope of the spatial power spectrum depends on both sonic and Alfven Mach numbers. This allows 
us to gauge the Alfven Mach number of the SMC and conclude that its gas pressure dominates over 
the magnetic pressure. The super- Alfvenic nature of the HI gas in the SMC is also highlighted by 
the bispectrum, a three-point correlation function which characterizes the level of non-Gaussianity 
in wave modes. We find that the bispectrum of the SMC HI column density displays similar large- 
scale correlations as numerical simulations, however it has localized enhancements of correlations. In 
addition, we find a break in correlations at a scale of ~ 160 pc. This may be caused by numerous 
expanding shells of a similar size. 

Subject headings: ISM: structure — MHD — turbulence: SMC 



1. INTRODUCTION 

In the recent decade, many advances in both observa- 
tions and computational models have provided new in- 
sights into the workings and evolution of the interstellar 
medium (ISM). The emerging picture is that interstel- 
lar turbulence plays the key role in ISM structure for- 
mation and evolution (McKee & Ostriker 2007). In the 
Galaxy and the Magellanic Clouds, the ISM is turbulent 
on scales ranging from less than a parsec to a few kilopar- 
secs (Crovisier & Dickey 1983; Stanimirovic et al. 1999; 
Deshpande et al. 2000; Dickey et al. 2001; Elmegreen 
et al. 2001; Elmegreen & Scalo 2004). Although the ob- 
servational evidence for the importance of turbulence in 
the ISM is overwhelming, many questions remain open. 
For example, what are the dominant energy sources and 
physical processes that convert kinetic energy into tur- 
bulence (Burkert 2006)? At what scales and through 
which modes is turbulent energy dissipated (Heyer & 
Zweibel 2004)? How do the level and type of turbulence 
depend on properties of the interstellar gas (e.g. pres- 
ence/absence of star formation, presence/absence of tidal 
effects, or the strength of magnetic field)? Since no com- 
plete theory of astrophysical turbulence exists, studying 
its effects on the multiphase ISM can be challenging and 
calls for a combination of numerical and observational 
efforts. 

Statistical studies have proved to be imp ortant in 
chara cterizing the magnetized turbulent ISM (jLazarianl 
2009), however the interpretation of results is not always 
straight forward. Several statistical methods have been 
extensively used for both observational and synthetic 
data. These statistics include probability density func- 
tions (PDFs), wavelets, the principal component analy- 
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sis, higher order moments, Velocity Coordinate Spectrum 
(VCS), and Velocity Channel Analysis (VCA), to name 
just a few (Gill & Henriksen 1990; Brunt & Heyer 2002; 
Kowal, Lazarian & Pogosyan 2000, 2004, 2006; Lazarian 
& Beresnyak 2007). 

Most of these statistical methods require large datasets 
with a large spatial or velocity dynamic range, and pro- 
duce a single, mostly one-dimensional, measure. This 
results in the lack of spatial information about turbulent 
properties across a given interstellar cloud, or a galaxy, 
making a connection with underlying physical properties 
highly difficult. 

In this paper we explore a new method for obtaining 
spatial information about the level and nature of ISM 
turbulence on the neutral hydrogen (HI) observations of 
the Small Magellanic Cloud (SMC). The SMC, a dwarf 
irregular galaxy in the Local Group, has a highly gas 
rich ISM environment (see, Stanimirovic et al. 1999, 
henceforth known as SX99), and is an exc ellent candi- 
date fo r ISM studies. Being nearby (60 kpc, Westerlund 
(1991)), the SMC is distant enough for all its objects to 
be treated as having roughly the same distance, unlike 
the Milky Way where distance determination is relatively 
uncertain. 

The HI observations of the SMC, obtained using the 
Australia Telescope Compact Array (ATCA) and the 
Parkes telescope (SX99), have been used for several in- 
vestigations, including the HI spatial power spectrum 
and the kinematic study of HI, which revealed the ex- 
istence of many expanding shells of gas and three su- 
pergiant shells. The power law index of HI density 
and velocity distributions wa s derived in SX99 and 
IStanimirovic &: Laz arian (2 0011 ). while the Genus statis- 
tic in lChepurnov et al.l (|2008f ) revealed spatial variations 
of HI morphology. Because the SX99 SMC data set is 
well studied, it is a perfect candidate to investigate new 
statistical methods. We can acquire new information, 
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but also test and confirm past results, as well as validate 
the promise of these statistical tools for further use in 
other observational studies. 

In this study, we investigate turbulent properties of 
the HI in the SMC by applying the higher order mo- 
ments on the HI column density image. We then use a 
database of MHD simulations to bootstrap the spatial 
distribution of the sonic Mach number across the SMC. 
The crucial aspect of our approach is the confluence of 
observations and numerical simulations: only by com- 
bining the two we can retrieve the spatial variations of 
turbulent properties. This is the reason why we oscillate 
between observational and synthetic data in this paper. 
We also investigate whether and how interstellar shocks 
leave footprints on the HI gas by employing the bispec- 
trum, a three point statistical measure, on the SMC HI 
column density image. Again, to interpret our results 
we apply the same statistics on the database of MHD 
simulated column density images. 

In particular, the paper is organized as follows. We 
start with § [2] by providing a brief summary of previ- 
ous work regarding the statistical methods used in our 
study. In § [3] we describe the SMC HI column density 
map and the database of numerical simulations of com- 
pressible MHD turbulence used for the comparison with 
observations. In § 4 we introduce higher order moments 
and their dependence on the sonic and Alfvenic Mach 
numbers. We then apply higher order moments on the 
SMC HI observations to derive an image of the sonic 
Mach number across the SMC, in § 5. We compare our 
results with an observational estimate of the sonic Mach 
number of the cold neutral medium (CNM) in the SMC, 
based on a comparison of the spin and kinetic tempera- 
ture of HI absorption profiles, in § |6] In § [7] we show how 
the power-law slope of the spatial power spectrum de- 
pends on the sonic Mach number and use this to gauge 
the Alfvenic Mach number of the HI gas seen in emis- 
sion. In §|S] we present an analysis of the bispectrum of 
the SMC, as well as a brief discussion of the noise and 
windowing effects. In § |9]we provide a discussion of our 
results, followed by our conclusions. 

2. BACKGROUND ON STATISTICAL METHODS USED IN 
THIS STUDY 

2.1. Higher order statistical moments 

Higher order statistical moments of density fluc- 
tuations have been studied extensively. For exam- 
ple, the variance of density fluctuations has been 
shown to increase with the sonic Mach num ber A4 S 
(|Nordlund & Padoanlfl999HOstriker et al.ll2001f ). There- 
fore, if turbulence is the dominant structuring mecha- 
nism, M s can be estimated from the variance of density 
fluctuations. However, the observable that is the most 
easily available from observations for various ISM trac- 
ers is the column density. While this is a less direct 
measure of turbulence compared to velocity, a compar- 
ison between observations and simulations is the most 
straight forward in case of column densities. 

Only very recently, Kowal et al. (2007), henceforth re- 
ferred to as KLB, have investigated how variance, skew- 
ness and kurtosis, the 2nd, 3rd and 4th order moments 
respectively, depend onM s . They found strong correla- 
tions: as the sonic Mach number increases, so does the 
Gaussian asymmetry of the column density (and density) 



PDFs due to gas compression via shocks, resulting in the 
increase of variance, skewness and kurtosis. KLB used 
limited resolution models of 128 3 in their study, while 
Burkhart et al. (2009), henceforth known as BFKL, saw 
the same trends using high resolution isothermal simu- 
lations. In both studies, the moments had little depen- 
dence on the Alfven Mach number (A4a), or the line- 
of-sight (LOS) orientation used for integrating 3-D sim- 
ulated data cubes. These studies are motivational, and 
imply that the sonic Mach number of turbulence in inter- 
stellar clouds could be characterized by variance, skew- 
ness and kurtosis of observed column density distribu- 
tions. 

2.2. Spatial power spectrum 

The two-dimensional spatial power spectrum charac- 
terizes the energy distribution over spatial scales. SX99 
found that the spatial power spectrum of individual HI 
velocity slices is well fit by a power-law, with an average 
slope of « —3. Stanimirovic & Lazarian (2001) showed 
that the power-law slope steepens when several channels 
are integrated together, and used this to estimate the 
density power-law slope of —3.3 and the velocity slope of 
—3.4. No obvious breaks or preferred scales were found 
over the range of 30 pc to 4 kpc. The density and velocity 
spectral slopes are similar, and in the case of incompress- 
ible MHD turbulence, density behaves as a passive scalar 
and thus scales in the same way as velocity (Monin & Ya- 
glom(1967), Lithwick & Goldreich (2001), Cho & Lazar- 
ian 2003). However, the estimated slopes are slightly 
more shallow than the predictions for the Kolmogorov 
spectrum ( fc~ n / 3 ), which is expected for incompressible 
fluids with a weak magnetic field. 

However, different types of turbulence are expected to 
have different spectral slopes. 

For example, in incompressible fluids with a strong 
magnetic field, the spectrum is expected to be even 
steeper and scale as fc~ 13 / 3 (Biskamp 2003). 

When the medium is supersonic (as we will see later ap- 
plies to parts of the SMC), these relations are no longer 
valid due to shocks forming highly asymmetric density 
structures. KLB demonstrated that the spectral slope 
of MHD turbulence becomes more shallow with increas- 
ing M. s . This can be understood as shocks in supersonic 
turbulence create more small scal e structure in density 
(|Beresnvak. Lazarian fc Choi l2005h . This behavior was 
found to be weakly dependent of the Alfven Mach num- 
ber. 

2.3. Bispectrum 

While the spatial power spectrum has long been ap- 
plied on both observations and simulations, it essen- 
tially uses only the amplitude of the Fourier transform 
of the initial signal, while the phase information is to- 
tally ignored. The bispectrum, however, is a three point 
statistical measure which incorporates both the ampli- 
tude and phase of the correlation of signals in Fourier 
space. Because of this, it can be used to search for 
nonlinear wave-wave interactions and characterize how 
shocks affect turbulent properties of the ISM. The bis- 
pectrum has been used in cosmology and gravitati onal 
wave studies to detect departures from Gau ssianity (|Frvl 
1998; Scoccimarroll2000l : iLiguori et al.ll2006l ). and for the 
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character ization of wave- wave interactions in labo ratory 
plasmas (jlntrator et all 119891 : I Tynan et all 1200 If ) . The 
first application of the bispectrum on synthetic astro- 
physical MHD turbulence was in BFKL. 

BFKL found a general correlation between the bispec- 
trum of 2D column density and 3D density maps for sim- 
ulated data cubes of 512 3 resolution. Also, supersonic 
models showed a much greater degree of correlation be- 
tween structures of different scales than subsonic mod- 
els. While comparing models with the same sonic Mach 
number, models with a stronger magnetic field (i.e. sub- 
Alfvenic) showed enhanced correlations. These results 
are encouraging and suggest that the bispectrum can be 
also used to provide insight into the nature of the turbu- 
lence cascade. 

3. DATA: SMC AND NUMERICAL SIMULATIONS 

3.1. HI Observations of the SMC 

The small-scale HI structure of the SMC was observed 
with the ATCA, a radio interferometer, in a mosaick- 
ing mode (Staveley-Smith et al. 1997). Observations of 
the same area were obtained also with the 64-m Parkes 
telescope, in order to map the distribution of large-scale 
features. Both sets of observations were then combined 
(see SX99), resulting in the final HI data cube with angu- 
lar resolution of 98", velocity resolution of 1.65 km s" 1 , 
and 1-er brightness temperature sensitivity of 1.3 K, to 
the continuous range of spatial scales between 30 pc and 
4.4 kpc. The velocity range covered with these observa- 
tions is 90-215 km s _1 . For details about the ATCA and 
Parkes observations, data processing, and data combina- 
tion (short spacings correction) see Staveley-Smith et al. 
(1997) and SX99. 

The HI column density image is shown in Figure [TJ 
We corrected the original image by multiplying the HI 
column density of each pixel (Njji in atom cm -2 ) with 
the correction factor f c derived by SX99: 



TABLE 1 
Description of the simulations 



MHD, 512 3 



fc 



l + 0.667(logiV ff /-21.4) 
1 



\ogN HI > 21.4 
\ogN HI < 21 A 



As the original Position-Position- Velocity (PPV) SMC 
cube has 578 x 610 x 78 pixels, the resulting column den- 
sity image is a 2D array with 578 x 610 pixels. 

3.2. A database of synthetic MHD cubes 

We generate a database of 13 3D numerical simula- 
tions of isothermal compressible (MHD) turbulence by 
using the code of KLB and varying the input values for 
the sonic and Alfvenic Mach number. The sonic Mach 
number is defined as A4 S = (|v|/C s ), where is v is the 
local velocity, C s is the sound speed, and the averaging is 
done over the whole simulation. Similarly, the Alfvenic 
Mach number is Ma = (W\/va), where va = |B|/^/p is 
the Alfvenic velocity, B is magnetic field and p is den- 
sity. KLB used resolution of 128 3 , while we use 512 3 . 
Details about KLB's code were published (see Cho et 
al. 2003) and the code was used in several studies. We 
briefly outline the major points of their numerical setup. 

The code is a second-order-acc urate hybrid essentially 
nonoscillatory (ENO) scheme (see lCho fc Lazarianl l2002f ) 
which solves the ideal MHD equations in a periodic box: 
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V x (v x B) = 0, (3) 



with zero-divergence condition V • B — 0, and an isother- 
mal equation of state p = C^p, where p is the gas pres- 
sure. On the right-hand side, the source term f is a 
random large-scale driving force 3 . We drive turbulence 
solenoidally, at wave scale k equal to about 2.5 (2.5 times 
smaller than the size of the box). This scale defines the 
injection scale in our models in Fourier space to minimize 
the influence of the forcing on the generation of density 
structures. Density fluctuations are generated later on by 
the interaction of MHD waves. The time t is in units of 
the large eddy turnover time (~ L/SV) and the length in 
units of L, the scale of energy injection. The scale of en- 
ergy dissipation is defined by the numerical diffusivity of 
the scheme 4 . The magnetic field consists of the uniform 
background field and a fluctuating field: B = B GXt + b. 
Initially b = 0. 

We divided our models into two groups corresponding 
to sub-Alfvenic (B cxt = 1-0) and super-Alfvenic (i? cx t — 
0.1) turbulence. For each group we computed several 
models with different values of gas pressure (see Table 
[I]). We ran compressible MHD turbulent models, with 
512 3 resolution, for t ~ 5 crossing times, to guarantee full 
development of energy cascade. Since the saturation level 
is similar for all models and we solve the isothermal MHD 
equations, the sonic Mach number is fully determined by 
the value of the isothermal sound speed, which is our 
control parameter. The models are listed and described 
in Table fU 

For each model, the results of MHD simulations are: 
the isothermal 3D density field, three components of ve- 
locity (V x , V y , V z ), and three components of magnetic 
field. As an example, Figure [2] (left) shows a density 
field for a sub-Alfvcnic, subsonic simulation. To calcu- 

3 f = pdv/dt 



vel y low 
s7j ri999l . 



— + V 

dt 



(pv)=0, (1) 



diffusion ones fsee lLiu Hi O shcr 1998; Levy, Puppo & R usso 
e.g.). The numerical diffusion depends not only on the adopted 
numerical scheme but also on the "smoothness" of the solution, 
so it changes locally in the system. In addition, it is also a time- 
varying quantity. 
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Fig. 1. — The final HI column density image after correction for self-absorption. The largest scale is ■ 
units of column density (cm -2 ). 



- 4.4 kpc. The intensity range is in 



late the column density distribution we integrate the 3D 
density fields along a given LOS. We introduce the fol- 
lowing nomenclature throughout the paper: "x column 
density" or "column density in the x direction" refers to 
the density cube being integrated along the x direction 
(parallel to the mean magnetic field). This description is 
similar for the y and z directions (perpendicular to B ext ). 
In the case of the Figure [2] (left), the LOS is along the z 
axis, and the magnetic field is oriented along the x axis. 

3.2.1. Scaling of the column density 

In order to compare simulations with the SMC data, we 
apply a simple scaling to the simulated column density 
distributions: 

N — (N) 



-^norrn ~ 



a(N) 



(4) 



where <j(N) denotes the standard deviation and (N) de- 
notes the arithmetic mean of the column density image. 
This method, often referred to as the standard score, 
standardizes all data sets used in this study. After the 
scaling, the new data set has values which represent the 
difference between the original data values and the sam- 
ple mean, in units of the standard deviation. The same 
scaling was applied on the SMC HI column density im- 
age. 

3.2.2. Simulating cloud boundaries 

Another consideration to take into account is the fact 
that simulated cubes are periodic and do not have bound- 
aries. In other words, simulations do not have a decrease 
in the column density values that is seen in the SMC 
data as one goes radially out from the center of the im- 
age. KLB showed that boundaries affect higher statis- 
tical moments of density, but have little effect on the 



spatial power spectrum and the bispectrum, generally 
only impacting the large scales. 

To introduce cloud boundaries to our simulations, we 
multiply the simulated 3D density fields by a 3D spher- 
ical function, which has a value of one within a sphere 
of radius R and is decaying outside of this radius with 
a Gaussian profile. We choose a R value of 128 pixels, 
which is 4 times smaller then the box size. An example 
of the simulation with a boundary is shown in the right- 
hand panel in Figure[2] Please note that we impose cloud 
boundaries on the simulated cubes of homogeneous tur- 
bulence after turbulence has been fully developed. This 
results in the simulated column density profiles being 
similar to the SMC profile, however, is different from 
simulations of turbulence being developed within cloud 
boundaries. 

4. HIGHER ORDER MOMENTS: SIMULATED DATA 

The most straightforward statistical properties of a dis- 
tribution are the mean value and the variance. The mean 
arithmetic value and the variance of the distribution of, 
for example, column density £ are given by: 



1 N 



and 



N 



1 N 



(5) 



(6) 



where N is the number of samples or points of the mesh 
in the case of simulation data. The mean value is a less 
important property in our studies so we do not consider 
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Fig. 2. — Left: An example of a simulated subsonic, sub-Alfvenic (M s = 0.7, M.A = 0.7) 512 3 density data cube demonstrating one 
possible direction of the LOS used to calculate the column density distribution. The mean magnetic field is oriented along the x direction. 
Here the LOS is along the z axis, although throughout the paper we integrate density fields along all three axes. Right: A supersonic 
column density image of a simulation with cloud boundaries, with a boundary radius similar (in pixel size) to what is seen for the SMC. 



it here; however, it is required to calculate higher mo- 
ments. Variance measures the width of the distribution 
and is always positby the third and fourth-order statis- 
tical moment. Skewness is defined as: 




If a distribution £ is Gaussian, its skewness is zero. Neg- 
ative skewness indicates the data are skewed in the left 
direction (the tail of the PDF is extended to the left, or 
towards low density values), while positive values imply 
that the distribution is skewed in the right direction (the 
tail of the PDF is extended to the right, or towards high 
density values). 

Kurtosis is a measure of whether a quantity has a dis- 
tribution that is peaked or flattened compared to a Gaus- 
sian distribution. Kurtosis is defined in a similar manner 
to skewness, only is derived from the forth order statis- 
tical moment. If a data set has positive kurtosis then 
the PDF of data values will have a distinct sharp peak 
near the mean and elongated tails. If a data set has a 
negative kurtosis then its PDF will be flat at the mean. 
Kurtosis is defined as: 

Variance, skewness and kurtosis of simulated column 
densities depend on turbulent properties, specifically the 
sonic Mach number. This is shown in Figure [3] for 
the simulated column density distributions with cloud 
boundaries. Error bars are computed by estimating the 
standard deviation between models with the same sonic 
Mach number but differing LOS orientations. All three 
moments depend almost linearly on the sonic Mach num- 
ber, for A4 S ;> 1 — 2. KLB noticed a similar trend for 
supersonic simulations but focused only on.M s = 0.2 — 2. 
Our simulations extend this result all the way to A4 S = 8. 
For Ai s <^ 0.5 KLB found that both skewness and kur- 
tosis have a rather flat dependence on the sonic Mach 
number, while the variance continues to decreases grad- 
ually for subsonic models. The increase of higher order 
moments with Ai s can be explained in that as the Mach 
number increases, the abundance of small-scale density 
structure increases, resulting in a broader density PDF. 



BFKL showed that the increase of higher-order moments 
with A4 S is not a resolution effect, since they replicated 
the result of KLB, (who used cubes of 128 3 resolution) 
with cubes of 512 3 resolution. 

Which one of three high-order moments represents the 
best statistics to describe turbulence in observational 
data? While variance has a linear dependence over a 
broad range of Ai s values, it depends on the exact scaling 
of the data set, making a direct comparison between sim- 
ulations and observations difficult. On the other hand, 
skewness and kurtosis are, by definition, normalized by 
the standard deviation and do not have this problem. 
In addition, variance changes with the length along the 
LOS, since perturbations can add up in a random walk 
fashion. All of this, as well as the result by KLB that 
kurtosis is the least affected by cloud boundaries, sug- 
gest that for our immediate comparison of observations 
with the MHD simulations, the higher order moments, 
skewness and kurtosis, are more appropriate statistics. 

4.1. Spatial Distribution of Skewness and Kurtosis 

While it is useful to know the Mach number of the ISM 
in a global sense, even more interesting is to see how it 
varies spatially and whether it correlates with local star 
forming regions where high turbulence is expected. To 
characterize small scale departures from Gaussian dis- 
tributions within the column density distributions, we 
create moment maps of the simulated column densities 
using a circular moving kernel. We do this by moving a 
circle of a given radius r across the image and calculating 
the skewness and kurtosis at all points. 

First we must decide on a kernel size that will en- 
close enough pixels to provide reliable statistics. Ac- 
cording to Tabachnick & Fidell (1996), the standard er- 
ror in skewness (SES) can be estimated by y/6/N and 
the standard error in kurtosis (SEK) can be estimated 

by y/24/N, where N is the number of points used to 
calculate the skewness/kurtosis. Generally, if a value of 
skewness/kurtosis falls within ±2x SES/SEK, then the 
distribution is considered to be normal, while values of 
skewness/kurtosis larger than twice the absolute value of 
SES/SEK imply significant non-Gaussianity. Figure [4] 
shows 2 x SES and SEK as a function of the kernel radius 
r. Clearly, SES/SEK is proportional to 1/r; the smaller 
the radius, the higher the absolute value of SES/SEK. 
Thus, we select r = 35 pixels, which corresponds to 
the point where 2x SES/SEK starts to deviate from 
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Fig. 3. — Variation of the variance (left), skewness (middle), and kurtosis (right) with the sonic Mach number derived from the final 
snapshot of each simulation with cloud boundaries applied. For all simulations Ma=0.7. Plotted data points are averages from data sets 
integrated along different LOS, with error bars representing the standard deviation between column density images with different LOS. As 
skewness and kurtosis increase almost linearly with M s , a linear function was fit over the range M s ~ 1 — 8 and is shown with a dashed-line. 
The measured skewness and kurtosis of the SMC HI column density image are: skewness= 1.52 ± 0.01 and kurtosis= 2.45 ± 0.02. Both 
of these values suggest that the HI in the SMC is supersonic, with a sonic Mach number > 3 within the error bars. This is discussed in 
Section 5. 



Skewness 
Kurtosis 



Log Radius of Beam 

Fig. 4. — The standard error in skewness (SES) and kurtosis (SEK) vs. radius of the circular kernel (related to the number of points 



while SEK= ^/2A/N 



We choose 



or pixels within the kernel used to calculate the skewness/kurtosis) . SES= y6/iV = \l 

r = 35 pixels to retain good resolution and low SEK/SES. These error measurement describes how large the statistics must be before they 
are deemed significant deviations from Gaussianity. If statistics are in this range then they are deemed generally Gaussian. If statistics are 
outside of this range then departures from Gaussianity have occurred. 



zero and ensures that Gaussianity is represented by 2x 
SES/SEK-0. 

Using a circular moving kernel we generate 3rd and 4th 
moment maps of simulated column density distributions. 
To ensure there are no resolution artifacts, we briefly ex- 
plore correlations between skewness and kurtosis of the 
derived moment maps on a pixel-to-pixel basis. Based on 
the work by KLB and BFKL, and also our Figure [3J we 
expect for supersonic models a correlation of both skew- 
ness and kurtosis with the sonic Mach number, and that 
skewness and kurtosis agree in sign and relative value. 
However, as subsonic models show unexpected behavior 
at low resolutions (from the KLB study), the dependence 
of skewness and kurtosis on the sonic Mach number is not 
strong. 

Figure [5] shows pixel-to-pixel comparison between 
skewness and kurtosis for two example models: a su- 
personic with M s = 7.0 (left panel) and a subsonic with 
Ai s = 0.7 (right panels). It is evident that the super- 
sonic model shows a good correlation between kurtosis 
and skewness, while the subsonic model does not. This 
is exactly as expected. This result proves that our ap- 
proach of deriving moment maps is not resolution limited 
although we use a smaller number statistics within each 
circular kernel (irr 2 — 1380 pixels). This also highlights 
another interesting property: regions of moment maps 



with a good correlation between skewness and kurto- 
sis could be interpreted by supersonic MHD turbulence, 
while a poor correlation may be caused by subsonic tur- 
bulence. We present in the appendix another technique 
that could be used to further identify subsonic regions. 



5. HIGHER ORDER MOMENTS: THE SMC HI COLUMN 
DENSITY IMAGE 

We first calculate the higher order moments of the 
whole SMC HI column density image. We find that the 
skewness is 1.52±0.01, and the kurtosis is 2.45 ± 0.02. 
The estimated uncertainties are determined by calculat- 
ing twice the SES or SEK. Figure [3] suggests that on 
average the HI in the SMC is supersonic with Ai s ^> 3. 

To characterize small scale departures from the Gaus- 
sian distribution within the SMC we create 3rd and 4th 
moment maps of the SMC using a circular kernel, as dis- 
cussed for simulations in Section H~T1 We repeat a simi- 
lar process for the SMC as we did with the simulations. 
We use a radius of r — 35 and calculate moments in a 
circular kernel moving across the SMC. This is equiva- 
lent to convolving the original HI image with a Gaussian 
function with a FWHM of 30'. The resultant moment 
images therefore have a lower resolution than the orig- 
inal HI image. After this we create a mask which cuts 
off pixels below a column density of 1.26 x 10 21 cm~ 2 in 
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Fig. 5. — Kurtosis vs. Skewness for isothermal simulations with a kernel of r = 35 pixels. The left is a supersonic model, while the right 
is a subsonic model. The supersonic model shows good correlation between skewness and kurtosis while the subsonic model does not. 



the HI column density image, to exclude increase in noise 
along the edges of the image. This is done to eliminate 
noisy pixels in moment maps, since pixels along the edge 
have significantly smaller column density values. 

Figures [6] shows the distribution of skewness and kur- 
tosis across the SMC (the HI column density image was 
standardized using the standard score method), with 
overlaid HI contours smoothed to 30' resolution. These 
maps retain the overall shape of the SMC, which can be 
seen in the HI contours. Familiar features of the SMC, 
such as the east wing and bar, regions of high star for- 
mation, can be picked out in these maps. 

Generally, skewness ranges from —0.5 to 1, with several 
discrete regions reaching higher positive or negative val- 
ues. For kurtosis, most pixels range from —1 to 1, with 
again several exceptions. Based on Figure 3, this suggest 
that large areas in the SMC have M s = — 2 and are 
subsonic or transonic. Along most of the bar the skew- 
ness and kurtosis correlate well. This suggests that the 
MHD turbulence could be the cause of local departures 
from Gaussianity in the HI column density distribution. 
We warn once again that both our study and earlier work 
find strong dependence of skewness and kurtosis on the 
sonic Mach number only for supersonic models, while 
the dependence for subsonic models is rather weak. In- 
terpreting subsonic turbulence is therefore difficult, and 
poor correlation between skewness and kurtosis could be 
due to either subsonic turbulence or some other physi- 
cal mechanisms. Several interesting regions stand out in 
Figure [H 



• Along the HI bar and the Eastern Wing region both 
skewness and kurtosis are negative, reaching values 
even lower than what is seen in Figure 3 for global 
averages. Based on expectations from numerical 
simulations this could be explained by the subsonic 
isothermal MHD turbulence, as shown in Figure 5. 
Most of the HI bar and the Wing therefore could be 
explained as being subsonic, with several discrete 
concentrations with almost no turbulence, proba- 
bly tracing quiescent regions and potential sites of 
future star formation. However, as subsonic regions 
are hard to constrain from simulations, additional 
processes may be also playing an important role. 



• From HI peaks radially outward, both skewness 
and kurtosis gradually increase. The highest values 
of both skewness (~ 2) and kurtosis (^3 — 4) are 
found along the HI bar and correspond to areas of 
compressed HI contours. These regions could be in- 
terpreted as having the highest level of turbulence 
(with M s =~ 4) in the SMC. This suggests that 
the most turbulent regions may be associated with 
the shearing flows and/or shocks between the bar 
and the surrounding HI. 

• Sudden change in the behavior of both skewness 
and kurtosis can be noticed around RA 01 h 10 m , 
Dec —72° 10'. For example, skewness flips from 
high positive to high negative values over an an- 
gular scale of ~ 60'. Interestingly, this flip hap- 
pens in the direction of an HI extension towards 
the LMC and again may be pointing out streaming 
or tidal motions caused by the interactions between 
the SMC and the LMC. 

In order to compare the range of values found in the 
SMC and the simulations we plot several histogram of the 
skewness and kurtosis values in Figure [7J These figures 
show that the majority of SMC pixels generally match 
well with the mildly supersonic models with M s ~ 1 — 
2. However, the SMC distribution appears to be more 
narrow in terms of both skewness and kurtosis values 
than the M s ~ 2 simulation. For subsonic and very 
supersonic models, the values of skewness and kurtosis 
of the SMC are not in the range of the simulations. 

5.1. The M s map of the SMC 

As shown above, with our limited angular resolution 
of 30' we see strong hints that different regions in the 
SMC have different turbulent properties, suggesting spa- 
tial variations oi M s - We would now like to combine 
all this information into an image of M s for the SMC. 
As already discussed, based on the KLB study, kurtosis 
is the least prone to effects caused by cloud boundaries 
and as Figure 3 suggests has an almost linear shape for 
M s ~ 1 — 8. Therefore, we fit a linear function for all 
values of kurtosis >^ —1. This function can be inverted 
to estimate M s : M s ~ (kurtosis + 1.44)/1.05. Pixels in 
the M s image with corresponding kurtosis < — 1 are set 
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Fig. 6. — (Left) Skewness of the HI column density derived using a circular kernel with r = 35 pixels. HI contours are overlaid. (Right) 
Kurtosis of the HI column density with the HI contours overlaid. The contours show the HI column density and range from 20 to 90%, 
with a step of 15%. 
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Fig. 7. — Histograms of the skewness and kurtosis maps from the simulations and the SMC HI column density image (see Figure [6), 
derived using the moving circle kernel with r = 35. The top row shows kurtosis, and the bottom row shows skewness. The SMC distributions 
are overplotted with dashed lines and agree best within the M s = 2.0 simulations. 



to zero and most likely mark subsonic regions, or could 
be caused by additional processes. The resultant map is 
shown in Figure [5] (left) . 

As expected already, most of the area of the SMC bar 
and the Eastern Wing could be interpreted as having sub- 
sonic or transonic Mach numbers, A4 S ^0 — 2. Several 
concentrated regions across the SMC indicate very qui- 
escent environments, most of them unfortunately have 



a size close to our angular resolution. Regions with the 
highest sonic Mach number, M s ~ 4, are found around 
the bar and correspond to compressed HI contours. 

We can quantify the fraction of HI with different 
A4 S . The most likely quiescent regions (set to zero in 
our map) comprise 8% of the mapped area. Regions 
with < Ai s < 1 comprise about 40%, while about 
the same fraction is contained in transonic regions with 
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Fig. 8. — (Left) The sonic Mach number image derived from the HI column density image of the SMC overlaid with the HI column 
density contours. The circle in the bottom-left shows the angular resolution of the image. (Right) Same as on the left but overlaid with 
the Ha contours smoothed to the same resolution. The Ha image is from Kennicutt et al. (1995). 



M. s = 1 — 2. Regions with higher turbulence, M s > 2 
constitute about 10% of the area. 

Figure [5] (right) shows contours of the Ha emission 
(obtained by Kennicutt et al. 1995) smoothed to 30' 
and overlaid on the A4 S map. Sites of the most recent 
star formation in the bar and the Eastern Wing have 
A4 S ~ 1, suggesting that the most turbulent regions are 
not associated with star formation. The most turbulent 
regions appear to trace shearing and tidal motions. We 
note though that our resolution of 30' is too low to trace 
individual star- forming regions. 

As an estimate of the velocity dispersion along the LOS 
is often used as a possible measure of HI turbulence, we 
investigated how the HI velocity dispersion compares to 
our A4 S map. Stanimirovic et al. (2004) found that re- 
gions with the highest dispersion (~ 30 km s _1 ) appear 
to be associated with the positions of the three largest 
super-giant shells. We find no obvious correlation be- 
tween the HI velocity dispersion and our A4 S map. 

This lack of correlation could be due to projection ef- 
fects, and/or large LOS depth of the SMC. In addition, 
velocity dispersion (similar to velocity centroids) is sub- 
ject to the interplay of several statistical effects (Lazar- 
ian & Esquivel 2003, Esquivel & Lazarian 2005, Esquivel 
et al. 2007). The most obvious is the contribution of 
both velocities and densities to the resulting measure. 
More subtle, but still important, are effects of velocity- 
density correlations and non-Gaussianity. These effects 
make centroids unreliable when dealing with supersonic 
turbulence (Esquivel et al. (2007). Effects of phase tran- 
sitions resulting in the existence of cold and warm HI and 
the pliable equation of state corresponding to interstellar 
HI are likely to act in a similar way. Thus the velocity 
dispersions obtained from spectral lines are substantially 
affected by densities, which are the quantities that we 
deal in this study. They should be distinguished from 
the true velocity dispersions, which, unfortunately, are 
not available through spectral line observations. 



6. COMPARISON WITH THE OBSERVATIONALLY 
INFERRED Ms OF THE CNM 

Constraining the sonic Mach number of the warm neu- 
tral medium (WNM) observationally is highly difficult 
due to the lack of direct measurements of the HI tem- 
perature. This is the main reason why in the previous 
section we used the HI column density image of the SMC 
and a database of numerical simulations to explore spa- 
tial variations of A4 S . To investigate reliability of our 
results we compare our derived A4 S distribution with 
results from an observational method which constrains 
temperature of the cold HI, and therefore Ai s of the 
CNM. This method is based on the comparison of the 
spin temperature (T s ) with the upper limit on the ki- 
netic temperature {Tk, m ax) for CNM clouds seen in HI 
absorption. 

For observational data we use HI absorption and emis- 
sion observations of 29 radio continuum sources be- 
hind the S MC by Dickey et al (2000). For each of 
the sources iDickev et alt (|2000l) provided the FWHM 
linewidth of the absorption spectra, and estimated the 
spin temperature of the CNM seen in absorption. By 
assuming Doppler broadening of HI velocity profiles, 
we can estimate the upper limit on the kinetic tem- 
perature as Tk. max = (FWHM/0.215) 2 . As shown in 
iHeiles fc Trolandl (|2003l ). the ratio of Tk^ max and the spin 
temperature, T s , is related to the ID mean square tur- 
bulent velocity: 

2 = m fTj^max _ A (g) 
m H \ T s J 

where k is the Boltzmann constant and uih is the mass 
of the hydrogen atom. Multiplying this by 3 gives the 
mean square 3D turbulence velocity V^ 3D of the CNM. 
To estimate the sonic Mach num ber we nee d to divide 
Vt,3D by the sound speed C s — \JkT s j fxrriH . We adopt 
a mean atomic weight of [i — 1.22 for the ISM of SMC 
abundance (see Mao et al. 2008). Using this, we can 
write: 
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For each of 29 radio continuum sources from 
iDickev et al.1 (|2000h we derive Tk, ma x and use the es- 
timated T s to calculate the sonic Mach number of the 
CNM clouds. IDickev et all (|2000ft used three different 
methods to estimate the spin temperature. As our aim 
is to compare M s with values derived from the LOS in- 
tegrated HI emission profiles, we use their LOS averaged 
spin temperature (or T ew )- The median value of the spin 
temperature is 43 K. 

A histogram of derived M s values is shown in Figure[5J 
The median Mach number for the whole sample of 29 
sources is M s = 4.7, and the histogram peaks around 
M s = 3.5 - 4. This suggests that the internal CNM 
macroscopic motions are highly supersonic. To compare 
these values with the ones derived using the method of 
higher-order statistical moments we show a histogram of 
data points from the M s map in the right panel of Fig- 
ure 03 While the observed histogram suffers from the 
low number statistics, due to a small number of suit- 
able (strong) sources, the two histograms have a similar 
shape: almost a Gaussian central distribution with a tail 
at higher Mach values. Obviously, the observationally in- 
ferred M s values suggest at least a factor of 2 higher level 
of turbulence across the SMC than what we derived using 
higher statistical moments. This is not surprising as the 
observed values trace predominately the CNM, while the 
M s map was derived using the HI column density im- 
age and therefore is affected by both CNM and WNM. 
We discuss this further in § \§\ In any case, to sample 
better CNM turbulence across the SMC future sensitive 
HI absorption observations of continuum sources behind 
the SMC are clearly needed. 

7. SPATIAL POWER SPECTRUM OF COLUMN DENSITY 

The slope of the Fourier transform of the two point 
autocorrelation function, or the spatial power spectrum, 
also provides information on important properties of tur- 
bulence flows. In particular, the slope of the power 
spectrum is known to depend on the sonic Mach num- 
ber. With a database of numerical simulations we ex- 
plore whether the slope of the spatial power spectrum 
depends also on the Alfven Mach number. The Alfven 
Mach number dependence would be especially interesting 
as we could use higher-order statistical moments to esti- 
mate the sonic Mach number, and from there the Alfven 
Mach number and the strength of the magnetic field. 

The power spectrum is defined as: 

P(k)= J2 A(k)-A*(k) (11) 

\k=k\ 

where k is the wavenumber and A(k) is the Fourier trans- 
form. Stanimirovic & Lazarian (2001) estimated the 
power-law slope of —3.3 for the spatial power spectrum 
of the HI column density image of the SMC. We measure 
the power spectrum slope for various simulated column 
densities. 

Figure [TU] shows how the power spectral slope changes 
with the sonic Mach number. The slope is increasingly 
shallow for supersonic models and levels off for very high 



Mach number turbulence. This is expected as higher 
Mach number turbulence has more density irregularities 
and more power on small scales. We have shown in this 
figure slopes for two values of Ma- 0.7 (dashed line) 
and 2.0 (dotted line). For both strong and weak mag- 
netic fields, the power spectrum slope increases with M Sl 
however this increase is steeper for super- Alfvenic and 
subsonic turbulence. For all cases the error bars were 
calculated from the standard deviation in the power spec- 
tral slope derived from column density images for three 
different LOS orientations. 

The SMC has a power spectrum slope of —3.3 and spa- 
tially mostly M s < 2. Based on the difference between 
slopes expected for sub- Alfvenic and super- Alfvenic tur- 
bulence in Figure [T7JJ it is likely that the super- Alfvenic 
description fits better the slope of the SMC. 

8. BISPECTRUM 

8.1. Calculation of the bispectrum 

The bispectrum is closely related to the power spec- 
trum. In a discrete system, the power spectrum is de- 
fined by Equation (11). In a similar way, the bispectrum 
is defined as: 

B(k 1 ,k 2 ) = J2 (12) 

ki—ki k~2— /C2 

where k\ and ki are the wave numbers of two interacting 
waves, and A(k) is the original discrete time series data 
with finite number of elements with A*{k) representing 

the complex conjugate of A(k). 

The bispectrum is a complex function which measures 
both phase and magnitude information between different 
wave modes. As this is the first application of the bispec- 
trum on the HI data we show in Figure [TT] a visual ex- 
ample of the difference between the power spectrum and 
the bispectrum. This figure shows the original HI column 
density image of the SMC (top left), and the same image 
but with manipulated phases (top right). To obtain the 
top right image we Fourier transformed the SMC image 
and randomized the phases with a Gaussian random dis- 
tribution but left the amplitudes the same. Even though 
the phases of the top images are very different, the power 
spectrum uses only amplitude information and is iden- 
tical for both images, as shown in the bottom panels of 
Figure 1111 However, the bispectrum looks very different 
for the two images, as expected, and offers an insight into 
the phase information. 

In practice, our calculation of the bispectrum involves 
performing a Fast Fourier Transform (FFT) of the col- 
umn density functions and the application of Equa- 
tion 1121 We randomly choose wavevectors and their di- 
rections, k\ and k 2 and iterate over them, calculating 
ka = k\ + k% . We limit the maximum length of the wave 
vectors to half of the box size. We normalize direction 
vectors to unity, calculate positions in Fourier space, and 
finally, compute the bispectrum, which yields a complex 
2D array. When plotting the bispectrum (Figures 12-15) 
we plot bispectral amplitudes, which give information 
about the degree of mode correlations in the original col- 
umn density distribution. 

8.2. Bispectrum of MHD Simulations 
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Fig. 9. — (Left) Histogram of the sonic Mach number of the CNM clouds in the SMC derived from the ratio Ta/T^^a^ using the HI 
absorption profiles against radio continuum sources. (Right) Histogram of the sonic Mach number from the Ms map derived using the 
higher statistical moments (see Figure [5J - 
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Fig. 10. — Power spectral slope vs. sonic Mach number for a sub-Alfvenic (Ma = 2.0, dotted line) and a super- Alfvenic (Ma = 0.7, 
dashed line) simulation. The spectral slope of the SMC (—3.3) is shown as a straight line. 



In order to interpret the bispectrum of the SMC HI 
column density we first compute the bispectrum of the 
MHD simulations (scaled by the standard score method 
and without applying cloud boundaries) and look for 
dependence of wave-wave correlations on the sonic and 
Alfven Mach numbers. Figure [T^] shows the bispectrum 
of simulated column density distributions for the case of 
fixed Ma = 0.7 and two extreme cases of M s : 0.7 and 
7.0. In Figure [13] we fix M s = 2.0 and look at two cases 
of M A - 0.7 and 2.0. 

The first thing to notice in the bispectral contour maps 
is that the k\ = k^ case always shows high correla- 
tion since this is a trivial case of two wave numbers 
being the same. For k\ ^ k^ the bispectral ampli- 
tude and isocontour shape vary with turbulent proper- 
ties but generally amplitude decreases gradually radially 
from k\ = &2 = 0. Models with circular ly/broad-shaped 
isocontours have high amplitudes and wave-wave corre- 
lations, while models with more narrow isocontours have 
lower amplitudes and therefore weaker correlations. For 
a fixed Ma, supersonic models show a higher degree of 
wave- wave correlations over subsonic models. For a fixed 
M s , models with a higher magnetic field (sub-Alfvenic, 
e.g. Ma =0.7) show somewhat stronger correlations 
then the models with a weaker magnetic field (super- 
Alfvenic) although this difference is not striking. 



8.3. Bispectrum of the SMC 

The bispectrum of the SMC HI column density is 
shown in Figurc [T4l The top row (first column) shows the 
bispectrum of the column density derived by integrating 
all 78 velocity channels. The middle and bottom rows 
show the bispectrum of the first (for channels 1 to 42) 
and the last half (for channels 43 to 78) of the data cube. 
To facilitate interpretation the wavenumbers ki and /c2 
are shown in terms of linear size 5 . While there are small 
differences in the bispectral amplitudes, there is essen- 
tially no significant difference in the contour shape for 
the three bispectra. 

To properly compare the SMC bispectrum with simula- 
tions we apply a windowing function on the SMC column 
density image to simulate the effect of periodic bound- 
aries. We use the Hanning window w(n) function, de- 
fined as: 

„(n) = 0^1-cos— j (13) 

where N is the number of pixels along the x or y axis, 
and n ranges from to N — 1. This windowing function 

5 The linear scale L = 4.4 kpc corresponds to the largest angular 
scale covered in the original image used to calculate the bispectrum, 
= 576 * 30" = 15,360". We scale then all wavenumbers k by 
2 X 15,360"/fc to express wavenumbers in terms of physical linear 
size at the distance of the SMC (60 kpc). 
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Fig. 11. — This figure illustrates differences between the bispectrum and the spatial power spectrum. The image at the top left is the 
original SMC HI column density image, while the image on the right has the same amplitudes however the phases have been randomized 
with a random Gaussian distribution. As expected, the power spectrum of the two images is identical, but the bispectrum shows a significant 
difference. 



M s =0J, M A =0.7 



Column 



ju=7.o, m,=oj 



Column D. -X Log(F 12 ) 





Fig. 12. — The amplitude of the bispectrum for the scaled simulated column density, derived by integrating along the x direction. The 
left plot shows a subsonic model, while the right plot is for a supersonic model. Both models have .Ma =0.7. These figures show the degree 
of correlation between wavenumbers fci and &2. The supersonic model has higher bispectral amplitudes, and more circular isocontours, 
therefore a stronger correlation between wave modes. 



makes the map periodic in the Fourier space. The second 
column of Figure [T4l shows the resultant bispectra for all 
three cases of HI column density integrations. The win- 
dowing seems to slightly change bispectral amplitudes 
on the large scales, but leaves the general structure of 
the bispectral contours intact, although some smoothing 
effects are observed. The increase in large scale correla- 
tions is at a level of ~ 10%. 

Before interpreting the SMC bispectrum we briefly in- 
vestigate how the presence of noise in the observational 
data affects the bispectrum. We made the noise image of 
the HI column density by taking a single velocity chan- 
nel from the SMC HI data cube without HI emission 
and integrating this image over the whole velocity range. 
The bispectrum of the noise image (normalized using the 



standard score method) is shown in Figure [T5] (right). 
In the case of purely Gaussian noise, we expect a random 
distribution over the whole space of wavelengths. As a 
result, the bispectrum will have a gradual increase in am- 
plitude, with the highest amplitude being at large k (or 
small scale). Figure [TS] (left) illustrates this effect for a 
pure Gaussian distribution. The bispectrum of the SMC 
noise image shows essentially the same trend, suggest- 
ing that we are dealing mainly with Gaussian (random) 
noise. As is obvious from Figures 12-14, the bispectral 
amplitude increases in the opposite direction (from large 
to small scales). The effect of noise is therefore the most 
important at the smallest spatial scales, and we keep this 
in mind when interpreting the bispectrum. 
Figure [JJ] shows the degree of wave- wave correlations 
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Fig. 13. — The amplitude of the bispectrum for the scaled simulated column density. The left plot shows a sub-Alfvenic model, while the 
right plot is for a super- Alfvenic model. Both models have J\4 S =2.0. These figures show the degree of correlation between wavenumbers 
k\ and k^. While the difference in amplitudes is less pronounced than in the previous figure, the model with a higher magnetic field (i.e. 
sub-Alfvenic) has more circular isocontours and therefore a slightly enhanced correlations between wave modes. 



in the HI distribution. Modes are obviously strongly cor- 
related at the largest scales and show weaker correlation 
at small scale. This trend is similar to what was found 
for various MHD simulations, and obviously very differ- 
ent from the case of Gaussian noise. In addition, the 
fci = &2 line shows the strongest correlations, as is also 
seen in the simulations. However, contrary to simula- 
tions where the level of correlation gradually decreases 
along this line, we see significant small-scale variations. 
Several local peaks are noticeable and a strong break at 
mid-scales. The break occurs around k\ — lt2 ~ 160 pc 
and is seen in all three integrations of the HI data cube. 
This could be due to a lack of interactions between tur- 
bulent eddies at the scale of ~ 160 pc, caused possibly by 
numerous expanding shells in the SMC. Hatzidimitriou 
et al. (2005) and Stanimirovic (2007) showed that shells 
in the SMC range in size from 30 to 800 kpc, however 
the radius distribution peaks at ~ 60 pc. The break in 
the bispectrum may be signifying the lack of correlations 
on scales close to the typical shell diameter, or some ad- 
ditional physical processes. 

We can also compare the bispectrum of the SMC with 
the bispectra of various simulations. By visual inspection 
of Figures 12-13, the closest simulation to the SMC in 
terms of contour shapes is the one for A4 S = 2.0 and 
M.A = 2.0. Further work in this area is required to derive 
a real measure to quantitatively compare observed and 
simulated bispectra. 

9. DISCUSSION AND SUMMARY 

9.1. Turbulence 'properties of the HI gas in the SMC 

Using 3rd and 4th statistical moments of the HI col- 
umn density image and boothstreping turbulent informa- 
tion from a database of isothermal MHD simulations, we 
have mapped spatial variations of the sonic Mach num- 
ber across the SMC. While most of the HI seen in emis- 
sion in the SMC appears to be subsonic or transonic, 
several supersonic regions have emerged from our study. 
It is interesting that these regions do not correlate well 
with the most recent sites of star formation and seem 
to point out to large scale shearing or tidal flows. Com- 
monly, it is believed that supernovae and superbubbles 
are the main drivers of galactic turbulence (McCray & 
Snow 1979), with a typical size of ~ 100 pc. While we do 
not have high enough resolution to see changes on such 
small scales (~ 10') in our derived map of M. s , most of 



the star-forming bar of the SMC appears to have sub- 
sonic or transonic properties when viewed at resolution 
of 30'. 

The most turbulent regions in the SMC may be tracing 
some kind of shearing flows between the SMC bar and the 
surrounding HI. This suggests that SMC's chaotic history 
with the LMC and our own Milky Way has probably left 
strong turbulent imprints on the HI gas. The lack of a 
turnover in the HI spatial power spectrum on the largest 
observed scales is also indicative of the fact that tur- 
bulent energy injection happe ns largely on scales larger 
than the size of the SMC (| Stanimirovic fc Lazarianl 
I2001fh Similarly, Goldman (2000) suggested that the 
HI turbulence in the SMC was induced by large-scale 
flows from tidal interactions with the Milky Way and 
the LMC about 2 x 10 s yrs ago. Such large-scale bulk 
flows could have generated turbulence through shear in- 
stabilities, and this turbulence has not have had enough 
time to decay. 

Most of the HI in the SMC has a sonic Mach number 
of 1-2. This is on average at least two times smaller 
than what we inferred from HI absorption observations 
for the CNM in the SMC, M s -3.5-4. Similarly, for the 
CNM in the Milky Way Heiles & Troland (2003) found 
M. s — 3 with a large dispersion. A sonic Mach number of 
about 4-5 is commonly assumed for cold gas in molecular 
clouds (Federrath et al. 2009). For example, Heyer et al. 
(2006) measured from CO observations M s = 4.2 ± 0.17 
for the Rosette molecular cloud, and Ai s = 4.7±0.12 for 
G216-2.5. On the other hand, Hill et al. (2008) found 
that the distribution of the warm ionized medium (WIM) 
in the Milky Way can be best fit by models for mildly 
supersonic turbulence with Ai s — 1.4 — 2.4. Turbulent 
properties of the HI in emission in the SMC are therefore 
closer to properties of the WIM in the Milky Way than 
properties of the CNM. This may suggest a large fraction 
of warm relative to cold HI being traced in HI emission. 

In the Milky Way, the HI is known to consist of at least 
two components with different temperature: the WNM 
with T warm = 6000 K and the CNM with T cold = 70 K. 
In addition, there could be a s ubstantial amount of ga s 
at intermediate temperatures (|Heiles fc Trolandl 12003). 
Due to its lower metallicity, the HI in the SMC has differ- 
ent properties. Dickey et al. (2000) found T co u — 40 K, 
in agreement with theoretical expectations by Wolfire 
et al. (1995) whereby the existence of the two-phase 
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Fig. 14. — The amplitude of the bispectrum for the SMC HI column density. The left column shows the bispectrum for the scaled 
unwindowed SMC image, while the right column shows the effects of a Hanning window on the bispectrum. The first row shows the SMC 
column density image derived by integrating over all 78 velocity channels. The second row is for the image derived by integrating over 
channels 1-42, while the third row shows the bispectrum of the image obtained by integrating the HI data cube over channels 43-78. The 
appropriate column density images used to derive the bispectrum are shown in the right-hand column. The bispectrum shows the degree 
of correlation between wavenumbers fci and ki in at varying depths of integration. We display the bispectral amplitude for spatial scales 
4400—80 pc as the bispectrum is very noisy at smaller scales. 



medium is possible only at higher pressures compared 
with the range that applies for solar neighborhood con- 
ditions. They also estimated the fraction of cold HI in 
the SMC to be <; 15%. This is lower than ~ 25% found 
for the Milky Way. 

As our simulations are isothermal it is obvious to won- 
der how does the multi-phase HI affect our statistics and 
conclusions. We investigate this in Figure [TBI by produc- 
ing a simulated data cube from a weighted combination 
of two cubes, one subsonic (.M s =0.7) and one supersonic. 
The subsonic cube represents contribution from warmer 
gas, while the supersonic cube represents the cold gas. 
We combined the two cubes with different emphasis on 
warm vs cold gas, obtained the column density image of 



the resultant cube, and calculated its moments. 

Figure [TBI shows that for the case when the supersonic 
cube has .M s =2.0 skewness of the final cube with up to 
50% of subsonic gas will be biased towards supersonic 
gas and will appear dominated by cold gas. If we in- 
crease the sonic Mach number of the supersonic cube 
to 4.0, the dominance of the higher turbulence is even 
more pronounced. A cube with up to 60-70% of subsonic 
gas and 25% of supersonic gas, will still have high skew- 
ness biased by the supersonic contribution. Considering 
that the HI column density image results in the mean 
A4 S = 1^2, significantly lower than what is expected 
for the cold HI, the fraction of the CNM along any LOS 
is most likely <; 25%. This supports the Dickey et al. 
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Fig. 15. — (Left) The amplitude of the bispectrum of a pure Gaussian distribution. It gradually increases from almost zero at small k to 
10 -3 at large k. (Right) The amplitude of the bispectrum derived for a single emission-less channel from the SMC HI data cube. 



(2000) estimate of the fraction of cold HI in the SMC 
being about 15%. 

9.2. Is the HI in the SMC sub-Alfvenic or 
Super- Alfvenic ? 

Two different statistical approaches in our study sug- 
gest that the HI gas in the SMC seen in emission is super- 
Alfvenic. As we have shown in Figure 1 101 in addition to 
the sonic Mach number the spectral slope of the spatial 
power spectrum is sensitive to the Alfvenic Mach number 
for Ai s < 2. The sub-Alfvenic models generally show 
steeper slopes due to large scale influence of magnetic 
fields. Thus, if one independently knows the sonic Mach 
number, it is possible to estimate the Alfvenic one us- 
ing just the column density data. While the dependence 
of the spectral slope on the Alfvenic Mach number has 
not received much attention in the past, it is somewhat 
expected. Essentially, magnetization decreases compres- 
sion in the shocks. Strong magnetic forces mix up density 
clumps preventing formation of isolated peaks, which re- 
sults in a steeper spectrum. In addition, in the sub- 
Alfvenic case we expect oblique shocks to be disrupted 
by Alfven sh earing, which in turn, produces m ore small 
scale shocks dBeresnvak. Lazarian fc Choll2005[) . 

Another indication that the HI in the SMC is super- 
Alfvenic comes from the bispectrum. The very sharp 
decrease in the bispectral amplitudes from large to small 
scales observed for the SMC is the closest to the trend 
found for simulated data for the case of A4 S ~ 2 and 
Ma ~ 2. Detailed comparison between simulated and 
observed bispectra awaits future work, however this qual- 
itative comparison is certainly encouraging. 

Assuming on average A4 S ~ 1 — 2, the power spec- 
trum slope suggests a super- Alfvenic HI in the SMC with 
M.A ~ 1 — 3. This is generally in agreement with the ob- 
servationally inferred strength of the magnetic field by 
Mao et al. (2008). Using their estimate for B ey ± — 2 /iG, 
a radius of the SMC of 2 kpc, the total hydrogen mass 
of 4.2 x 10 9 M Q , and a typical velocity dispersion of 20 
km s _1 (Stanimirovic et al. 2004), we estimate Ma ~ 3. 
As the Alfvenic Mach number shows the nature of the 
interplay between gas pressure and magnetic fields, it ap- 
pears that the gas pressure in the SMC dominates over 
the magnetic pressure. 

9.3. Intermittency in mode correlations? 



Our bispectrum analysis of the SMC HI data was the 
first attempt to apply bispectrum on observed astrophys- 
ical data. While more detailed comparison between ob- 
servations and simulations awaits future work, we clearly 
see trends in the bispectral amplitudes similar to what 
was found for simulations of supersonic MHD turbulence. 
The most interesting finding is, however, the effect of 
small-scale variations in the fei = hi correlations and a 
strong break in correlations at a scale of 160 pc. Such 
small-scale variations, or jumps, have not been seen in 
the bispectrum of simulated data. We can speculate 
about several possible scenarios that could explain their 
existence. The jumps could be caused by the energy in- 
jection due to processes other than turbulence affecting 
specific spatial scales. Alternatively, the jumps may be 
marking the presence of colder or multi-phase gas. Simi- 
larly, the observed break in the bispectrum at about 160 
pc is intriguing. As we already pointed out, it is interest- 
ing that most expanding shells in the SMC (more than 
500 were cataloged so far) have a diameter of ~ 120 pc. 
The break could be due to the lack of correlations on 
scales similar to the distance between two shell centers. 
Obviously this will require further studies. 

9.4. Limitations of the present study 

A natural question to ask is how results presented in 
this paper depend on the resolution of numerical sim- 
ulations. For example, Kritsuk et al. (2007) investi- 
gated how resolution of numerical simulations affects the 
power spectrum of density. These authors found that 
the spectral index estimates based on low resolution sim- 
ulations bear large uncertainties due to the bottle neck 
contamination, and that the power spectra of 512 3 simu- 
lations are substantially shallower then models with res- 
olution of 1024 3 . However, while Kritsuk et al. (2007) 
only examined hydrodynamic turbulence, Beresnyak et 
al. (2008) showed that the slopes were very different be- 
tween the MHD and pure hydrodynamic cases. For in- 
stance, the slopes for hydrodynamic simulations showed 
a pronounced and well defined bottleneck effect, while 
the MHD slopes were much less affected. This is indica- 
tive of MHD turbulence being less local than the Kol- 
mogorov turbulence, and suggests that our simulations 
will be less affected by resolution. In addition, Kritsuk 
et al. (2007) found a difference in the slope between hy- 
drodynamic 512 3 and 1024 3 simulations to be 0.17. This 
would result in a change of dM s ~ 0.5 only and will not 
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Fig. 16. — Skewncss vs. percent of gas that is subsonic for column density of a 512 3 cube weighted with supersonic and subsonic gas. 
Supersonic gas generally dominates the skewness of the column density PDF. 



change our interpretation. We also add that in the case 
of higher statistical moments and the bispectrum BFKL 
confirmed trends noticed by KLB at lower resolution of 
128 3 . 

Another issue that should be further addressed and 
that could affect our results is the type of numerical 
forcing of turbulence. Federrath et al. (2009) recently 
investigated the effects of solenoidal vs. compressive 
(divergence-free vs. curl-free) forcing on a variety of 
statistics including PDFs and higher order moments. 
They found that both types of driving mechanisms are 
compatible with observations of molecular clouds how- 
ever, depending on the data studied, one type could be 
superior then the other in terms of the statistics and re- 
produced observables. This implies that different regions 
in the SMC may exhibit statistical signatures of cither 
compressive or solenoidally driven turbulence. 

9.5. Summary 

We have investigated a new method for constraining 
turbulent properties of the ISM, specifically the sonic 
Mach number, by using the HI column density image and 
a database of numerical simulations with a range of sonic 
and Alfvenic Mach numbers. By applying the 3rd and 
4th statistical moments on both observed and simulated 
data we have derived the spatial distribution of the sonic 
Mach number across the SMC with angular resolution of 
30'. To provide an estimate of the Alfvenic Mach number 
we used two approaches: the spatial power spectrum and 
the bispectrum. Using the database of numerical simu- 
lations we have shown that the spatial power spectrum 
varies with both the sonic and Alfvenic Mach numbers. 
If the sonic number is known the Alfvenic number can 
be constrained from this dependence. The bispectrum 
shows the level of correlation between turbulent eddies 
of different size and depends greatly on the sonic Mach 
number, and somewhat on the Alfvenic Mach number. 
By comparing the bispectra of observations and simula- 
tions we have gauged the importance of magnetic fields 
relative to the gas pressure in the SMC. The following 
results were discussed in the paper. 



• Skewness and kurtosis of the HI column density 
generally correlate well and are within the range ex- 
pected from MHD simulations. This suggests that 
departures from Gaussianity could be interpreted 
as being governed by MHD turbulence. 

• Most of the HI in the SMC bar and the Eastern 
Wing is subsonic or transonic with M s ~ — 2. 
Sites of most recent star formation have M s ~ 1- 
Regions with the highest skewness and kurtosis, 
which could be interpreted as having M s ~ 4, cor- 
respond to the edges of the bar. The most turbu- 
lent regions are most likely tracing tidal or shearing 
flows. The fraction of the SMC with different tur- 
bulent properties is: 10% with M s > 2, 80% with 
< M s < 2, and about 10% with very low values 
of M a . 

• Using HI absorption profiles from Dickey et al. 
(2000) we have estimated that the CNM in the 
SMC is highly supersonic with M s = 3.5 — 4. This 
is at least a factor of two higher than what we mea- 
sured from the higher statistical moments for the 
HI gas seen in emission. One possible reason for 
this discrepancy could be that HI emission is dom- 
inated by warm gas and the fraction of the CNM 
in the SMC is £ 20%. 

• The slope of the spatial power spectrum and the 
bispectrum suggest that the HI in the SMC is 
super- Alfvenic with Ma ~ 1 — 3. This is implies 
that the gas pressure dominates over the magnetic 
pressure. 

• The bispectrum of the HI column density shows 
large scale wave correlations suggesting a large 
scale energy injection mechanism. Contrary to sim- 
ulations which show a smooth decrease of wave- 
wave correlations from large to small scales, the 
SMC bispectrum shows localized enhancements 
of correlations and at least one prominent break 
at <~ 160 pc. We speculate that the multi- 
phase medium, and/or energy injection by pro- 
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cesses other than turbulence, could be responsible 
for the correlation jumps. The break on the other 
hand appears at a scale similar to the diameter of 
the majority of expanding shells in the SMC. 
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APPENDIX 

Here we present a technique that could be used to further illuminate subsonic regions in the ISM. 
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Fig. 17. — Kurtosis vs. Skewness for isothermal simulations. The left is subsonic while the right is supersonic. These figure are similar 
to [5] except that the beam was modified with 7 extra points set to zero. The supersonic models remain relatively unchanged in their trend 
while subsonic models with the modified distribution obtain a very tight anti-correlation due to extra points pushing them away from 
Gaussianity. The supersonic model is unaffected because it is already very positively skewed. This could be employed as an additional 
method to detect subsonic gas using moments 
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Fig. 18. — Kurtosis and Skewness map with beams modified with 7 zero points side by side. Areas where they do not agree in sign (dark 
blue in skewness and dark blue/red in kurtosis) point to anti-correlation , which indicates regions that are potentially subsonic. These 
regions correspond to areas of low kurtosis and skewness in the previous maps, yet are more clearly seen here. The largest of these subsonic 
regions is seen along the bridge between the bar and the east wing (See Figure [5}. 

In Figure [T7] we plot kurtosis vs. skewness in a manner similar to Figure [5j The difference between this figure 
and Figure [5] is that we added seven zero points into the beam (instead of 3841 points there are now 3848 points). 
This produces almost no change in the supersonic models since they already have high positive skewness and kurtosis. 
Supersonic models still show strong correlation between skewness and kurtosis. 

However a big change is seen in the subsonic skewness vs. kurtosis. The additional zero points shift the distribution 
from Gaussian (the mean of our simulations with no scaling applied is unity) to negative skewness and very peaked 
kurtosis producing a tight anti-correlation. Because this technique only strongly affects subsonic areas we can use it 
to locate subsonic turbulence in the SMC with skewness and kurtosis by looking for anti-correlation and very high 
values. These properties also hold for simulations with cloud boundaries imposed. Note that one must use caution 
here and carefully examine the distribution of data. 
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The application of this technique to the SMC data is shown in Figure [TH This plot shows the skewness and 
kurtosis maps of the SMC side by side with the modified beam. Indeed, most regions are unchanged from the analysis 
in Section [5l However a few regions stick out with signatures that are subsonic, that is, anti-correlation between 
skewness and kurtosis. Two regions of the highest kurtosis are located in the area between the HI bar and the Eastern 
Wing. While kurtosis reaches values of 4-8, the corresponding skewness values are negative, ~ — 1. Again, such 
combination of skewness and kurtosis may correspond to simulations of subsonic isothermal MHD turbulence, and/or 
points out additional processes at work. 
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